Simulations with a single SNP & different sample sizes

1 Setup

First, simulate genotype data. Then, for each subject’s genotype, simulate replicate traits.

Code

library(magrittr)

Code

set.seed(2023-04-15)
n_tr <- 10000
n_te <- 100
n <- n_tr + n_te
# n is total number of subjects
n_reps <- 50
# n_reps is number of replicates
n_snp <- 1
# n_snp is number of SNPs
# simulate genotypes matrix
geno <- matrix(sample(c(0,1,2), n*n_snp, replace=TRUE, prob=c(0.25,0.5,0.25)), nrow=n, ncol=n_snp)
# simulate phenotype
beta <- 1
y <- as.numeric(geno %*% beta) %*% t(rep(1, n_reps)) + 
    matrix(data = rnorm(n * n_reps), nrow = n, ncol = n_reps)
# prepare for splitting into training and testing
# get test subject ids
test_ids <- sample(1:n, n_te, replace = FALSE)
# organize data
dat <- tibble::as_tibble(y) %>%
    dplyr::rename_with(function(x){ num <- stringr::str_extract(x, "[0-9]+")
                                    return(paste0("pheno", num))}
                        ) %>%
    dplyr::bind_cols(geno %>% tibble::as_tibble() %>% dplyr::rename(geno = 1)) %>%
    dplyr::mutate(id = 1:n) %>% # fix this when using more than one SNP
    dplyr::mutate(in_test_set = id %in% test_ids)

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.

Code

# split into training and testing
training <- dat %>% dplyr::filter(!in_test_set)
testing <- dat %>% dplyr::filter(in_test_set)
testing2 <- testing %>% dplyr::mutate(fold = as.integer(NA))
# use all training with leave K out
alpha <- 0.1

Above, we specified the number of replicates for the simulations. We created 50 replicate traits for the same 1.01^{4} subjects. Note that each subject has only 1 SNP(s).

2 python code, adapted from Barber’s code for JK+ paper

Code

import numpy as np
import pandas as pd

def leastsq_minL2(X,Y,X1,tol=1e-8):
    uX,dX,vX = np.linalg.svd(X)
    rX = (dX>=dX[0]*tol).sum()
    betahat = (vX[:rX].T/dX[:rX]).dot(uX[:,:rX].T.dot(Y))
    return X1.dot(betahat)

def compute_PIs(X,Y,X1,alpha = 0.1, K=10, fit_muh_fun = leastsq_minL2):
    n = len(Y)
    #n1 = X1.shape[0]
    n1 = len(X1)
    ###############################
    # CV+
    ###############################
    n_K = np.floor(n/K).astype(int)
    base_inds_to_delete = np.arange(n_K).astype(int)
    resids_LKO = np.zeros(n)
    muh_LKO_vals_testpoint = np.zeros((n,n1))
    for i in range(K):
        inds_to_delete = (base_inds_to_delete + n_K*i).astype(int)
        X0 = np.delete(X, inds_to_delete)
        X0_2d = X0[:, np.newaxis]
        Y0 = np.delete(Y, inds_to_delete)
        X0d = X[inds_to_delete]
        X0d_2d = X0d[:, np.newaxis]
        X1_2d = X1[:, np.newaxis]
        muh_vals_LKO = fit_muh_fun(X0_2d, Y0, np.r_[X0d_2d, X1_2d])
        resids_LKO[inds_to_delete] = np.abs(Y[inds_to_delete] - muh_vals_LKO[:n_K, 0])
        for inner_K in range(n_K):
            muh_LKO_vals_testpoint[inds_to_delete[inner_K]] = muh_vals_LKO[n_K:, 0]
    ind_Kq = (np.ceil((1-alpha)*(n+1))).astype(int)
    PIs_dict = {'CV+' : pd.DataFrame(\
                    np.c_[np.sort(muh_LKO_vals_testpoint.T - resids_LKO,axis=1).T[-ind_Kq], \
                        np.sort(muh_LKO_vals_testpoint.T + resids_LKO,axis=1).T[ind_Kq-1]],\
                           columns = ['lower','upper'])}
    return pd.concat(PIs_dict.values(), axis=1, keys=PIs_dict.keys())

Code

# simulation
# simulation
n_vals = np.array([100, 1000, 10000, 100000]) # training set size
n1 = 100 # test set size
SNR = 10
ntrial = 50
alpha = 0.1
#dim_vals = np.arange(5,205,5)
d = 1
beta = np.random.normal(size=d)
beta = beta/np.sqrt((beta**2).sum()) * np.sqrt(SNR)
X1 = np.random.normal(size=(n1,d))
Y1 = X1.dot(beta) + np.random.normal(size=n1)
method = 'CV+'
# define new objects for storing Y and X
Y_dic = dict.fromkeys(n_vals, 0)
X_dic = dict.fromkeys(n_vals, 0)
# loop over n_vals
for n in n_vals:
    results = pd.DataFrame(columns = ['itrial','d','method','coverage','width'])
    Xa = np.zeros((n, ntrial))
    Ya = np.zeros((n, ntrial))
    for itrial in range(ntrial):
        X = np.random.normal(size=(n,d))
        Xa[:, itrial] = X[:,0]
        Y = X.dot(beta) + np.random.normal(size=n)
        Ya[:, itrial] = Y
        # store X, Y for later use with R code
        # store as numpy arrays
        PIs = compute_PIs(X = X, Y = Y, X1 = X1)
        coverage = ((PIs[method]['lower'] <= Y1)&(PIs[method]['upper'] >= Y1)).mean()
        width = (PIs[method]['upper'] - PIs[method]['lower']).mean()
        results.loc[len(results)]=[itrial,d,method,coverage,width]
    X_dic[n] = Xa
    Y_dic[n] = Ya

Code

saveRDS(py$PIs, file = "PIs.rds")

Code

res_list <- list()
vars_list_list_list <- list()
n_folds_vec <- c(5, 10, 20, 50)
for (n_folds in n_folds_vec){
    # partition training data into K folds
    folds <- split(training$id, sample(rep(1:n_folds, length.out = n_tr)))
    training2_pre <- training %>% 
        dplyr::mutate(fold = id %>% purrr::map_int(~which(sapply(folds, function(x) . %in% x)))) %>%
        dplyr::arrange(fold)
    tictoc::tic() # timing
    tl <- list()
    n_per_fold_vec <- c(n_tr / n_folds, n_tr / (n_folds * 10))
    vars_list_list <- list()
    for (n_per_fold in n_per_fold_vec){
        training2 <- training2_pre %>%
            dplyr::group_by(fold) %>%
            dplyr::slice_sample(n = n_per_fold) %>%
            dplyr::ungroup()
        # store each trait's outputs
        out <- list()
        vars_list <- list()
        # loop over traits
        for (trait_num in 1:n_reps){
            tr2_one_trait <- training2 %>%
                dplyr::select(id, fold, geno, tidyselect::ends_with(paste0("pheno", trait_num))) %>%
                dplyr::rename(pheno = tidyselect::ends_with(paste0("pheno", trait_num)))
            te2_one_trait <- testing2 %>%
                dplyr::select(id, fold, geno, tidyselect::ends_with(paste0("pheno", trait_num))) %>%
                dplyr::rename(pheno = tidyselect::ends_with(paste0("pheno", trait_num)))
            
            # regress leaving one fold out
            preds <- list()
            vars <- list()
            for (fold_num in 1:n_folds) {
                # get training data
                train <- tr2_one_trait %>% dplyr::filter(fold != fold_num)
                # get testing data
                test <- tr2_one_trait %>% dplyr::filter(fold == fold_num)
                # fit model
                fit <- lm(pheno ~ 1 + geno, data = train)
                # predict
                foo <- test %>% dplyr::bind_rows(te2_one_trait)
                foo$pred <- predict(fit, newdata = foo)
                foo$fold_left_out <- fold_num
                result <- foo %>%
                    dplyr::mutate(beta1_hat = coef(fit)[2],
                                beta0_hat = coef(fit)[1],
                                se_beta1_hat = summary(fit)$coefficients[2, 2],
                                se_beta0_hat = summary(fit)$coefficients[1, 2]
                    )
                # save predictions
                preds[[fold_num]] <- result
                te_geno_mat <- cbind(1, te2_one_trait$geno)
                tr_geno_mat <- cbind(1, train$geno)
                vars[[fold_num]] <- diag(te_geno_mat %*% solve(t(tr_geno_mat) %*% tr_geno_mat) %*% t(te_geno_mat))
            }
            vars_list[[trait_num]] <- vars
            # assemble predicted values
            # get absolute residuals
            preds_training <- preds %>%
                dplyr::bind_rows() %>%
                dplyr::filter(!is.na(fold)) %>% # keep only training data
                dplyr::mutate(absolute_residual = abs(pheno - pred)) %>%
                dplyr::select( - fold_left_out)
            preds_test <- preds %>%
                dplyr::bind_rows() %>%
                dplyr::filter(is.na(fold))
            # get indexes
            plus_index <- ceiling((1 - alpha) * (nrow(preds_training) + 1))
            minus_index <- floor(alpha * (nrow(preds_training) + 1))
        
            # go one by one through test set (testing2)
            test_list <- list()
            for (i in 1:nrow(testing2)){
                tt <- testing2[i, ]
                pt2 <- preds_test %>% 
                    dplyr::filter(id == tt$id) %>% # our only use of tt
                    dplyr::rename_with(function(x)paste0("test_", x)) 
                    # pt2 contains the five predicted values for a single test subject
                nrow(pt2) # 5
                preds_all <- preds_training %>%
                    dplyr::left_join(pt2, by = c("fold" = "test_fold_left_out")) %>%
                    dplyr::mutate(test_fitted_plus_absolute_residual = test_pred + absolute_residual, 
                                test_fitted_minus_absolute_residual = test_pred - absolute_residual) 
                uu <- sort(preds_all$test_fitted_plus_absolute_residual)[plus_index]
                ll <- sort(preds_all$test_fitted_minus_absolute_residual)[minus_index]
                # make a tibble with exactly one row
                test_list[[i]] <- preds_all %>%
                    dplyr::select(test_id, test_geno, test_pheno, test_beta1_hat, fold) %>%
                    dplyr::mutate(lower = ll, upper = uu) %>%
                    dplyr::distinct() %>%
                    tidyr::pivot_wider(names_from = fold, 
                                        values_from = test_beta1_hat,
                                        names_prefix = "beta1_hat_fold_"
                                        )
            }
            test_tib <- test_list %>%
                dplyr::bind_rows() %>%
                dplyr::mutate(in_interval = test_pheno >= lower & test_pheno <= upper) %>%
                dplyr::mutate(interval_width = upper - lower) %>%
                dplyr::mutate(training_set_size = n_per_fold * n_folds,
                                trait_num = trait_num)
            out[[trait_num]] <- test_tib
        }
        tl[[as.character(n_per_fold * n_folds)]]  <- out
        vars_list_list[[as.character(n_per_fold * n_folds)]] <- vars_list
    }
    res_list[[as.character(n_folds)]] <- tl
    vars_list_list_list[[as.character(n_folds)]] <- vars_list_list
    tictoc::toc() # timing
}

189.042 sec elapsed
193.453 sec elapsed
200.736 sec elapsed
224.502 sec elapsed

Code

saveRDS(res_list, "res_list.rds")
saveRDS(vars_list_list_list, "vars_list_list_list.rds")

3 Organize results

Code

#test_tib_thin <- test_tib %>%
#    dplyr::select(test_id, test_geno)
tt_intermediate <- tl %>%
    dplyr::bind_rows(.id = "id") 
results <- tt_intermediate %>%
    dplyr::group_by(training_set_size, trait_num) %>%
    dplyr::summarise(mean_interval_width = mean(interval_width),
                    sd_interval_width = sd(interval_width),
                    mean_in_interval = mean(in_interval),
                    sd_in_interval = sd(in_interval), 
                    beta1_hat_fold_1 = mean(beta1_hat_fold_1),
                    beta1_hat_fold_2 = mean(beta1_hat_fold_2),
                    beta1_hat_fold_3 = mean(beta1_hat_fold_3),
                    beta1_hat_fold_4 = mean(beta1_hat_fold_4),
                    beta1_hat_fold_5 = mean(beta1_hat_fold_5),
                    median_interval_width = median(interval_width)
                    ) %>%
                    dplyr::ungroup() %>%
                    dplyr::mutate(mean_b1 = purrr::pmap_dbl(.l = list(beta1_hat_fold_1,
                                                                        beta1_hat_fold_2, 
                                                                        beta1_hat_fold_3, 
                                                                        beta1_hat_fold_4,
                                                                         beta1_hat_fold_5), 
                                                           .f = function(x, y, z, w, v) mean(c(x, y, z, w, v))),
                                    sd_b1 = purrr::pmap_dbl(.l = list(beta1_hat_fold_1,
                                                                        beta1_hat_fold_2, 
                                                                        beta1_hat_fold_3, 
                                                                        beta1_hat_fold_4,
                                                                         beta1_hat_fold_5), 
                                                           .f = function(x, y, z, w, v) sd(c(x, y, z, w, v)))
                    ) 

`summarise()` has grouped output by 'training_set_size'. You can override using
the `.groups` argument.

Code

results %>%
    knitr::kable() %>%
    print()

training_set_size	trait_num	mean_interval_width	sd_interval_width	mean_in_interval	sd_in_interval	beta1_hat_fold_1	beta1_hat_fold_2	beta1_hat_fold_3	beta1_hat_fold_4	beta1_hat_fold_5	median_interval_width	mean_b1	sd_b1
1000	1	3.269663	0.0020488	0.90	0.3015113	0.9065156	0.9121913	0.9060416	0.9048478	0.9188715	3.268853	0.9096936	0.0058607
1000	2	3.344751	0.0023533	0.91	0.2876235	1.0445959	1.0464540	1.0443695	1.0332784	1.0392876	3.343471	1.0415971	0.0053573
1000	3	3.274158	0.0010671	0.94	0.2386833	0.9813570	0.9765217	0.9883275	0.9702371	0.9831070	3.274411	0.9799101	0.0068596
1000	4	3.230431	0.0016232	0.86	0.3487351	0.9111773	0.9119977	0.9151600	0.9135225	0.9184691	3.229089	0.9140653	0.0028941
1000	5	3.286248	0.0010091	0.92	0.2726599	0.9771232	0.9775981	0.9772230	0.9771884	0.9765326	3.287091	0.9771331	0.0003834
1000	6	3.320046	0.0015877	0.91	0.2876235	0.9946676	0.9978031	1.0119408	0.9960080	1.0008850	3.318530	1.0002609	0.0069322
1000	7	3.296412	0.0000000	0.89	0.3144660	0.9582067	0.9469858	0.9592551	0.9581995	0.9497622	3.296412	0.9544819	0.0056777
1000	8	3.317076	0.0016047	0.91	0.2876235	1.0228708	1.0077899	0.9995220	0.9989276	1.0127770	3.318388	1.0083775	0.0099684
1000	9	3.382647	0.0035541	0.90	0.3015113	0.9246682	0.9303803	0.9239508	0.9242849	0.9226916	3.384131	0.9251952	0.0029920
1000	10	3.290401	0.0007914	0.85	0.3588703	0.9839049	0.9978354	0.9718226	0.9894680	0.9905806	3.290028	0.9867223	0.0096936
1000	11	3.406227	0.0020223	0.92	0.2726599	1.0551749	1.0710386	1.0711884	1.0673222	1.0589285	3.406510	1.0647305	0.0073005
1000	12	3.279119	0.0032105	0.93	0.2564324	1.1133192	1.1120771	1.1019529	1.1167999	1.0992799	3.278216	1.1086858	0.0076258
1000	13	3.279855	0.0004555	0.84	0.3684529	1.0272876	1.0214734	1.0331090	1.0402648	1.0277830	3.280171	1.0299836	0.0070712
1000	14	3.221249	0.0007762	0.85	0.3588703	0.9793727	0.9827520	0.9806434	0.9683182	0.9839642	3.221994	0.9790101	0.0062382
1000	15	3.304753	0.0024335	0.91	0.2876235	1.0406638	1.0280019	1.0373577	1.0310384	1.0272879	3.303581	1.0328699	0.0058970
1000	16	3.351328	0.0016736	0.90	0.3015113	1.0630697	1.0605694	1.0561969	1.0552951	1.0522724	3.351343	1.0574807	0.0043107
1000	17	3.339365	0.0075022	0.92	0.2726599	1.0488649	1.0595641	1.0578711	1.0498131	1.0581815	3.344983	1.0548589	0.0050902
1000	18	3.388968	0.0055995	0.89	0.3144660	0.9920304	0.9821338	0.9725179	0.9795348	0.9800220	3.389309	0.9812478	0.0070302
1000	19	3.234607	0.0019073	0.91	0.2876235	1.0390109	1.0389052	1.0323969	1.0266685	1.0271897	3.236107	1.0328343	0.0060223
1000	20	3.281298	0.0006824	0.89	0.3144660	0.9269577	0.9225596	0.9040631	0.9020235	0.9139623	3.280792	0.9139132	0.0109922
1000	21	3.249487	0.0003837	0.91	0.2876235	0.9732818	0.9754162	0.9797883	0.9740106	0.9758393	3.249844	0.9756672	0.0025254
1000	22	3.243039	0.0012789	0.84	0.3684529	1.0176006	1.0150738	1.0252889	1.0219341	1.0264490	3.244161	1.0212693	0.0048804
1000	23	3.338817	0.0036488	0.93	0.2564324	0.9772634	0.9640830	0.9729947	0.9631500	0.9686649	3.341840	0.9692312	0.0059684
1000	24	3.425091	0.0037725	0.91	0.2876235	1.0430569	1.0346992	1.0301639	1.0403048	1.0304766	3.425764	1.0357403	0.0057934
1000	25	3.282211	0.0063019	0.94	0.2386833	1.0470885	1.0346177	1.0349041	1.0429263	1.0443080	3.276165	1.0407689	0.0056866
1000	26	3.213914	0.0017614	0.92	0.2726599	0.9783352	0.9903133	0.9872459	0.9844983	0.9827820	3.213518	0.9846349	0.0045332
1000	27	3.361788	0.0028798	0.86	0.3487351	0.9303066	0.8989762	0.9201872	0.9251130	0.9236144	3.360140	0.9196395	0.0121121
1000	28	3.379681	0.0023469	0.90	0.3015113	0.9623537	0.9656788	0.9583561	0.9563788	0.9583019	3.380762	0.9602139	0.0037503
1000	29	3.377388	0.0022226	0.95	0.2190429	1.0268036	1.0402019	1.0375707	1.0408578	1.0284271	3.376673	1.0347722	0.0066728
1000	30	3.342649	0.0035026	0.91	0.2876235	0.9871109	0.9863959	0.9989223	0.9974864	0.9843121	3.339364	0.9908455	0.0068148
1000	31	3.293695	0.0004773	0.83	0.3775252	0.9277105	0.9324002	0.9427483	0.9241295	0.9257188	3.293649	0.9305415	0.0074974
1000	32	3.333125	0.0047400	0.90	0.3015113	0.9819737	0.9903822	0.9906274	1.0054269	0.9960257	3.335702	0.9928872	0.0086253
1000	33	3.324350	0.0021101	0.93	0.2564324	1.0466477	1.0511812	1.0496322	1.0506528	1.0393048	3.322694	1.0474837	0.0048971
1000	34	3.302055	0.0021719	0.91	0.2876235	0.9796539	0.9925029	0.9815801	0.9880037	0.9858494	3.302621	0.9855180	0.0051228
1000	35	3.312941	0.0043681	0.93	0.2564324	1.0056364	1.0019997	1.0080133	1.0110617	0.9930155	3.313399	1.0039453	0.0069516
1000	36	3.317381	0.0024905	0.91	0.2876235	0.9635208	0.9649013	0.9629495	0.9610914	0.9588560	3.315525	0.9622638	0.0023449
1000	37	3.321484	0.0028642	0.91	0.2876235	1.0256105	1.0296355	1.0360205	1.0385688	1.0457209	3.323085	1.0351112	0.0078357
1000	38	3.291529	0.0023839	0.87	0.3379977	1.0279565	1.0339062	1.0369824	1.0349016	1.0320923	3.292975	1.0331678	0.0034065
1000	39	3.283259	0.0014321	0.93	0.2564324	1.0222746	1.0293725	1.0190563	1.0315141	1.0233214	3.283279	1.0251078	0.0051737
1000	40	3.254464	0.0013634	0.88	0.3265986	1.0535078	1.0519398	1.0618976	1.0659103	1.0481263	3.254882	1.0562764	0.0073706
1000	41	3.310142	0.0020953	0.93	0.2564324	1.0414822	1.0385111	1.0437353	1.0400287	1.0442206	3.311724	1.0415956	0.0024213
1000	42	3.296851	0.0026634	0.87	0.3379977	0.9895380	0.9838824	0.9892810	0.9767553	0.9979014	3.299058	0.9874716	0.0078110
1000	43	3.285878	0.0006534	0.87	0.3379977	0.9526910	0.9505892	0.9448929	0.9606618	0.9462730	3.286320	0.9510216	0.0062453
1000	44	3.395559	0.0037983	0.93	0.2564324	0.9509533	0.9426644	0.9480655	0.9336279	0.9281187	3.397338	0.9406859	0.0096377
1000	45	3.384050	0.0009418	0.90	0.3015113	1.0281364	1.0346181	1.0401679	1.0469562	1.0417212	3.384337	1.0383199	0.0071934
1000	46	3.361692	0.0012464	0.84	0.3684529	1.0500323	1.0631663	1.0451493	1.0532041	1.0730478	3.362904	1.0569199	0.0111683
1000	47	3.385481	0.0039291	0.87	0.3379977	0.9316182	0.9322810	0.9346191	0.9302921	0.9173239	3.386566	0.9292268	0.0068363
1000	48	3.324511	0.0021389	0.88	0.3265986	0.9783534	0.9858048	0.9770990	0.9870862	0.9822730	3.326377	0.9821233	0.0044063
1000	49	3.263428	0.0031921	0.88	0.3265986	1.0090399	1.0239101	1.0106192	1.0135401	1.0133401	3.265303	1.0140899	0.0058062
1000	50	3.383038	0.0030075	0.92	0.2726599	0.9909738	1.0001257	0.9973402	0.9833335	0.9952424	3.384251	0.9934031	0.0065455
10000	1	3.270154	0.0004563	0.89	0.3144660	0.9702581	0.9696415	0.9725096	0.9666922	0.9664444	3.269968	0.9691092	0.0025549
10000	2	3.307498	0.0006781	0.90	0.3015113	1.0068604	1.0082770	1.0059084	0.9994799	1.0064373	3.307266	1.0053926	0.0034202
10000	3	3.305666	0.0004471	0.94	0.2386833	0.9952225	0.9907088	0.9930630	0.9954499	0.9957652	3.305828	0.9940419	0.0021456
10000	4	3.297982	0.0004889	0.88	0.3265986	0.9905126	0.9920318	0.9917018	0.9938894	0.9901000	3.298420	0.9916471	0.0014882
10000	5	3.311855	0.0002613	0.92	0.2726599	0.9896307	0.9899452	0.9920244	0.9902287	0.9901152	3.312098	0.9903888	0.0009417
10000	6	3.299189	0.0001528	0.91	0.2876235	0.9866049	0.9861679	0.9884676	0.9901826	0.9866258	3.299099	0.9876098	0.0016891
10000	7	3.294630	0.0009200	0.89	0.3144660	1.0159411	1.0133711	1.0149591	1.0154214	1.0162966	3.294148	1.0151979	0.0011406
10000	8	3.290669	0.0002296	0.90	0.3015113	1.0051638	1.0038946	1.0043200	1.0072521	1.0082875	3.290570	1.0057836	0.0019053
10000	9	3.288675	0.0009101	0.87	0.3379977	0.9800655	0.9792149	0.9779191	0.9764331	0.9820672	3.289365	0.9791400	0.0021359
10000	10	3.294845	0.0006079	0.85	0.3588703	0.9928709	0.9917809	0.9889710	0.9926706	0.9953541	3.295227	0.9923295	0.0022988
10000	11	3.344140	0.0004833	0.92	0.2726599	0.9951009	0.9952136	0.9929867	0.9919202	0.9957589	3.344597	0.9941961	0.0016537
10000	12	3.268679	0.0002491	0.93	0.2564324	1.0210084	1.0186432	1.0193218	1.0222108	1.0224976	3.268835	1.0207364	0.0017126
10000	13	3.301849	0.0001953	0.85	0.3588703	1.0024792	1.0048416	1.0054082	1.0095662	1.0053171	3.302035	1.0055224	0.0025560
10000	14	3.277762	0.0003858	0.86	0.3487351	1.0145859	1.0152086	1.0133260	1.0109958	1.0160803	3.278137	1.0140393	0.0019751
10000	15	3.250141	0.0005377	0.90	0.3015113	0.9795953	0.9815929	0.9785730	0.9821062	0.9805872	3.250432	0.9804909	0.0014417
10000	16	3.249044	0.0004771	0.89	0.3144660	0.9985478	0.9973796	0.9991926	0.9985771	0.9990079	3.249487	0.9985410	0.0007058
10000	17	3.291970	0.0004443	0.91	0.2876235	1.0105684	1.0107225	1.0080642	1.0088449	1.0113155	3.292373	1.0099031	0.0013794
10000	18	3.314226	0.0003154	0.88	0.3265986	0.9867850	0.9857631	0.9865365	0.9867038	0.9881312	3.313972	0.9867839	0.0008551
10000	19	3.278342	0.0007176	0.91	0.2876235	1.0007617	1.0007086	1.0018955	1.0020897	1.0013254	3.277835	1.0013562	0.0006330
10000	20	3.296240	0.0002258	0.89	0.3144660	0.9693706	0.9709751	0.9676739	0.9677549	0.9748950	3.296422	0.9701339	0.0029862
10000	21	3.274367	0.0001712	0.92	0.2726599	1.0084111	1.0099171	1.0057049	1.0077447	1.0086067	3.274410	1.0080769	0.0015424
10000	22	3.268773	0.0007709	0.85	0.3588703	0.9778421	0.9767741	0.9783067	0.9803107	0.9781879	3.268425	0.9782843	0.0012838
10000	23	3.312616	0.0007980	0.92	0.2726599	1.0180404	1.0196089	1.0170376	1.0190702	1.0164472	3.313095	1.0180408	0.0013289
10000	24	3.257415	0.0002974	0.90	0.3015113	0.9955714	0.9952546	0.9984949	0.9976871	0.9959567	3.257691	0.9965929	0.0014190
10000	25	3.295426	0.0002744	0.94	0.2386833	1.0181587	1.0156793	1.0150684	1.0111794	1.0142163	3.295356	1.0148604	0.0025277
10000	26	3.299758	0.0005659	0.93	0.2564324	1.0004877	0.9988031	1.0017621	0.9971864	1.0004475	3.300080	0.9997374	0.0017711
10000	27	3.319314	0.0018348	0.85	0.3588703	0.9938729	0.9890072	0.9927564	0.9911029	0.9888362	3.318052	0.9911151	0.0022325
10000	28	3.246731	0.0005051	0.90	0.3015113	0.9897918	0.9934104	0.9913962	0.9922393	0.9920222	3.246702	0.9917720	0.0013255
10000	29	3.307375	0.0002463	0.96	0.1969464	0.9833423	0.9846720	0.9834209	0.9847303	0.9846258	3.307402	0.9841583	0.0007105
10000	30	3.261481	0.0003443	0.90	0.3015113	0.9897589	0.9928072	0.9883256	0.9902908	0.9905425	3.261813	0.9903450	0.0016222
10000	31	3.290697	0.0002159	0.82	0.3861229	0.9922776	0.9958630	0.9971857	0.9948750	0.9956903	3.290675	0.9951783	0.0018212
10000	32	3.244219	0.0003519	0.90	0.3015113	0.9980821	1.0018115	0.9987914	0.9999459	1.0005863	3.244546	0.9998435	0.0014701
10000	33	3.305079	0.0002525	0.92	0.2726599	0.9957653	0.9933296	0.9917444	0.9944223	0.9959041	3.305276	0.9942331	0.0017454
10000	34	3.304851	0.0005965	0.91	0.2876235	0.9895069	0.9962579	0.9943861	0.9943288	0.9932770	3.305276	0.9935513	0.0025030
10000	35	3.281931	0.0005446	0.93	0.2564324	0.9684046	0.9681545	0.9675809	0.9725219	0.9681103	3.281415	0.9689544	0.0020167
10000	36	3.291591	0.0008047	0.92	0.2726599	1.0045561	1.0022102	1.0023347	1.0041388	1.0032860	3.290827	1.0033052	0.0010489
10000	37	3.286436	0.0003196	0.91	0.2876235	1.0096608	1.0070383	1.0104226	1.0109595	1.0099677	3.286549	1.0096098	0.0015187
10000	38	3.292834	0.0001629	0.87	0.3379977	0.9992385	1.0001884	0.9967046	1.0015196	0.9996198	3.292883	0.9994542	0.0017635
10000	39	3.275772	0.0005378	0.93	0.2564324	1.0072665	1.0123215	1.0116986	1.0097169	1.0096925	3.276295	1.0101392	0.0019890
10000	40	3.267455	0.0001610	0.89	0.3144660	1.0026875	0.9972174	1.0022171	1.0019432	1.0002842	3.267612	1.0008699	0.0022335
10000	41	3.275336	0.0006660	0.92	0.2726599	0.9793571	0.9795859	0.9790515	0.9740876	0.9764112	3.274790	0.9776987	0.0023896
10000	42	3.307056	0.0003436	0.86	0.3487351	1.0095064	1.0082879	1.0087154	1.0058256	1.0090916	3.307382	1.0082854	0.0014471
10000	43	3.285535	0.0003932	0.86	0.3487351	1.0164067	1.0162693	1.0148107	1.0183307	1.0158239	3.285173	1.0163283	0.0012822
10000	44	3.307926	0.0003645	0.92	0.2726599	1.0042092	0.9999930	0.9992875	1.0017402	0.9998453	3.307590	1.0010150	0.0020079
10000	45	3.334554	0.0007744	0.90	0.3015113	0.9996387	0.9985265	1.0006972	0.9968591	0.9978923	3.334912	0.9987228	0.0014949
10000	46	3.317126	0.0004997	0.84	0.3684529	1.0054035	1.0014818	0.9999203	1.0028222	1.0091165	3.316666	1.0037488	0.0036127
10000	47	3.267832	0.0002567	0.86	0.3487351	0.9959679	0.9906880	0.9923052	0.9910029	0.9933180	3.267895	0.9926564	0.0021293
10000	48	3.298945	0.0002295	0.87	0.3379977	0.9994844	1.0003358	1.0000943	0.9975878	0.9985393	3.298980	0.9992083	0.0011408
10000	49	3.270811	0.0006287	0.88	0.3265986	0.9828106	0.9834019	0.9791488	0.9849327	0.9813058	3.271317	0.9823200	0.0021978
10000	50	3.295095	0.0003907	0.92	0.2726599	0.9744410	0.9798955	0.9773824	0.9741856	0.9790768	3.295288	0.9769962	0.0026131

4 Figures

4.1 Boxplots for interval width

Code

library(ggplot2)
tt_intermediate %>%
    ggplot(aes(y = interval_width, colour = as.factor(training_set_size), x = as.factor(trait_num)))  +
    geom_boxplot()

Code

#ggsave(here::here("figures", "interval_width_boxplot.png"), width = 10, height = 10)

4.2 Relationship between \(\hat\beta\) and median interval width

Code

p1 <- results %>%
    ggplot(aes(x = mean_b1, y = median_interval_width, colour = as.factor(training_set_size), replicate_num = trait_num)) +
    geom_point()
plotly::ggplotly(p1, tooltip = c("x", "y", "colour", "replicate_num"))

5 Session Info

Code

sessioninfo::session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.0 (2023-04-21)
 os       Ubuntu 18.04.6 LTS
 system   x86_64, linux-gnu
 ui       X11
 language en_US:
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Detroit
 date     2023-05-07
 pandoc   1.19.2.4 @ /usr/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 cli           3.6.1   2023-03-23 [1] CRAN (R 4.2.3)
 colorspace    2.1-0   2023-01-23 [1] CRAN (R 4.2.2)
 crosstalk     1.2.0   2021-11-04 [1] CRAN (R 4.1.2)
 data.table    1.14.8  2023-02-17 [1] CRAN (R 4.2.2)
 digest        0.6.31  2022-12-11 [1] CRAN (R 4.2.2)
 dplyr         1.1.1   2023-03-22 [1] CRAN (R 4.2.3)
 ellipsis      0.3.2   2021-04-29 [2] CRAN (R 4.2.1)
 evaluate      0.20    2023-01-17 [1] CRAN (R 4.2.2)
 fansi         1.0.4   2023-01-22 [1] CRAN (R 4.2.2)
 farver        2.1.1   2022-07-06 [1] CRAN (R 4.2.3)
 fastmap       1.1.1   2023-02-24 [1] CRAN (R 4.2.3)
 generics      0.1.3   2022-07-05 [1] CRAN (R 4.2.3)
 ggplot2     * 3.4.2   2023-04-03 [1] CRAN (R 4.2.3)
 glue          1.6.2   2022-02-24 [1] CRAN (R 4.2.0)
 gtable        0.3.3   2023-03-21 [1] CRAN (R 4.2.3)
 here          1.0.1   2020-12-13 [2] CRAN (R 4.1.1)
 htmltools     0.5.5   2023-03-23 [1] CRAN (R 4.2.3)
 htmlwidgets   1.6.2   2023-03-17 [1] CRAN (R 4.2.2)
 httr          1.4.5   2023-02-24 [1] CRAN (R 4.2.3)
 jsonlite      1.8.4   2022-12-06 [1] CRAN (R 4.2.3)
 knitr         1.42    2023-01-25 [1] CRAN (R 4.2.3)
 labeling      0.4.2   2020-10-20 [2] CRAN (R 4.0.3)
 lattice       0.21-8  2023-04-05 [1] CRAN (R 4.2.3)
 lazyeval      0.2.2   2019-03-15 [2] CRAN (R 4.0.3)
 lifecycle     1.0.3   2022-10-07 [1] CRAN (R 4.2.2)
 magrittr    * 2.0.3   2022-03-30 [1] CRAN (R 4.2.0)
 Matrix        1.5-4   2023-04-04 [1] CRAN (R 4.2.3)
 munsell       0.5.0   2018-06-12 [2] CRAN (R 4.0.3)
 pillar        1.9.0   2023-03-22 [1] CRAN (R 4.2.3)
 pkgconfig     2.0.3   2019-09-22 [2] CRAN (R 4.0.3)
 plotly        4.10.1  2022-11-07 [1] CRAN (R 4.2.2)
 png           0.1-8   2022-11-29 [1] CRAN (R 4.2.3)
 purrr         1.0.1   2023-01-10 [1] CRAN (R 4.2.2)
 R6            2.5.1   2021-08-19 [2] CRAN (R 4.1.1)
 rappdirs      0.3.3   2021-01-31 [2] CRAN (R 4.0.3)
 Rcpp          1.0.10  2023-01-22 [1] CRAN (R 4.2.2)
 reticulate    1.28    2023-01-27 [1] CRAN (R 4.2.3)
 rlang         1.1.0   2023-03-14 [1] CRAN (R 4.2.2)
 rmarkdown     2.21    2023-03-26 [1] CRAN (R 4.2.3)
 rprojroot     2.0.3   2022-04-02 [2] CRAN (R 4.2.0)
 scales        1.2.1   2022-08-20 [1] CRAN (R 4.2.3)
 sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.1.2)
 stringi       1.7.12  2023-01-11 [1] CRAN (R 4.2.2)
 stringr       1.5.0   2022-12-02 [1] CRAN (R 4.2.3)
 tibble        3.2.1   2023-03-20 [1] CRAN (R 4.2.3)
 tictoc        1.1     2022-09-03 [1] CRAN (R 4.2.1)
 tidyr         1.3.0   2023-01-24 [1] CRAN (R 4.2.3)
 tidyselect    1.2.0   2022-10-10 [1] CRAN (R 4.2.2)
 utf8          1.2.3   2023-01-31 [1] CRAN (R 4.2.3)
 vctrs         0.6.1   2023-03-22 [1] CRAN (R 4.2.3)
 viridisLite   0.4.1   2022-08-22 [1] CRAN (R 4.2.3)
 withr         2.5.0   2022-03-03 [1] CRAN (R 4.2.0)
 xfun          0.38    2023-03-24 [1] CRAN (R 4.2.3)
 yaml          2.3.7   2023-01-23 [1] CRAN (R 4.2.3)

 [1] /net/mulan/home/fredboe/R/x86_64-pc-linux-gnu-library/4.0
 [2] /net/mario/cluster/lib/R/site-library-bionic-40
 [3] /usr/local/lib/R/site-library
 [4] /usr/lib/R/site-library
 [5] /usr/lib/R/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /net/mulan/home/fredboe/miniconda3/envs/barber1/bin/python3
 libpython:      /net/mulan/home/fredboe/miniconda3/envs/barber1/lib/libpython3.10.so
 pythonhome:     /net/mulan/home/fredboe/miniconda3/envs/barber1:/net/mulan/home/fredboe/miniconda3/envs/barber1
 version:        3.10.9 (main, Mar  8 2023, 10:47:38) [GCC 11.2.0]
 numpy:          /net/mulan/home/fredboe/miniconda3/envs/barber1/lib/python3.10/site-packages/numpy
 numpy_version:  1.23.5
 
 python versions found: 
  /net/mulan/home/fredboe/miniconda3/envs/barber1/bin/python3
  /net/mulan/home/fredboe/miniconda3/envs/barber1/bin/python
  /net/mulan/home/fredboe/miniconda3/bin/python
  /net/mulan/home/fredboe/miniconda3/envs/ldsc/bin/python

──────────────────────────────────────────────────────────────────────────────

Code

# git commit info
gr <- git2r::repository(here::here()) %>%
    git2r::commits()
gr[[1]] 

[af93d90] 2023-05-07: fix: fixed typo